##########################################################################
# HANGARTNER, LAUDERDALE, SPIRIG (2025) INFERRING INDIVIDUAL PREFERENCES #
##########################################################################

#####################################################
########### BAYESIAN SINGLE YEAR ANALYSIS ###########
#####################################################################
# Define function to estimate preferences given fixed decision rule #
#####################################################################

# Y is a vector of M binary decisions, with 1 corresponding to decision in favor of the asylum seeker
# Ij is a M x 3 factor matrix with the three judges (columns), all of which have the same factor levels

####################
## mixture models ##
####################

stancode_covariate_mixture = '

data {
int<lower=1> M;	
int<lower=1> N;
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
real<lower=0,upper=1> delta;
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector<lower=0,upper=1>[2] pimat[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M){	// loop over cases

// Define chair submodel
pimat[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);

// Define median submodel
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);
thetamats3[j] = rep_vector(0.0,3);        
pimat[j,2] = 0.0;

// if all judges are known
thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);
thetamats3[j] = sort_asc(thetamat3[j]);
pimat[j,2] = thetamats3[j,2];


// Define mixture model from submodels
pi[j] = pimat[j,1]*delta + pimat[j,2]*(1-delta);

} // end loop over cases


}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

vector[3] thetaalt;	
vector[3] thetaalts;
real pialt;

real inconsistj[M];	
real inconsistjj[M];	
real inconsist;	
real pt1;
real pt2;

real pialt_median;
real inconsistj_median[M];	
real inconsistjj_median[M];	
real inconsist_median;	

real pialt_chair;
real inconsistj_chair[M];	
real inconsistjj_chair[M];	
real inconsist_chair;	

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);

// simulate probability of different outcome under alternative panels
        for (jj in 1:M){ 
            thetaalt[1] = inv_logit(logittheta[Ij[jj,1]] + phi[Lj[j]]);
            thetaalt[2] = inv_logit(logittheta[Ij[jj,2]] + phi[Lj[j]]);
            thetaalt[3] = inv_logit(logittheta[Ij[jj,3]] + phi[Lj[j]]);
            thetaalts = sort_asc(thetaalt);
            pialt = thetaalt[1]*delta + thetaalts[2]*(1-delta); // define predicted probability from sub-models
            pialt_chair = thetaalt[1];
            pialt_median = thetaalts[2];
            pt1 = (Y[j]*(pi[j] > pialt)*((pi[j] - pialt)/(pi[j])));
            pt2 = ((1-Y[j])*(pi[j] < pialt)*((pialt - pi[j])/(1 - pi[j])));
            inconsistjj[jj] = (pt1 + pt2);
            
            pt1 = (Y[j]*(pi[j] > pialt_chair)*((pi[j] - pialt_chair)/(pi[j])));
            pt2 = ((1-Y[j])*(pi[j] < pialt_chair)*((pialt_chair - pi[j])/(1 - pi[j])));
            inconsistjj_chair[jj] = (pt1 + pt2);
            
            pt1 = (Y[j]*(pi[j] > pialt_median)*((pi[j] - pialt_median)/(pi[j])));
            pt2 = ((1-Y[j])*(pi[j] < pialt_median)*((pialt_median - pi[j])/(1 - pi[j])));
            inconsistjj_median[jj] = (pt1 + pt2);
        }
    inconsistj[j] = mean(inconsistjj);
    inconsistj_chair[j] = mean(inconsistjj_chair);
    inconsistj_median[j] = mean(inconsistjj_median);
}
ll = sum(llj);
inconsist = mean(inconsistj);
inconsist_chair = mean(inconsistj_chair);
inconsist_median = mean(inconsistj_median);


}

'



stancode_covariate_mixture_split = '

data {
int<lower=1> M;	
int<lower=1> N;
int<lower=0,upper=1> SC[M];
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
real<lower=0,upper=1> delta1;
real<lower=0,upper=1> delta0;
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector<lower=0,upper=1>[2] pimat[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M){	// loop over cases

// Define chair submodel
pimat[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);

// Define median submodel
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);
thetamats3[j] = rep_vector(0.0,3);        
pimat[j,2] = 0.0;

// if all judges are known
thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);
thetamats3[j] = sort_asc(thetamat3[j]);
pimat[j,2] = thetamats3[j,2];

// Define mixture model from submodels, split by covariate
if (SC[j]==1) { // if cov = 1
pi[j] = pimat[j,1]*delta1 + pimat[j,2]*(1-delta1);
} else {
pi[j] = pimat[j,1]*delta0 + pimat[j,2]*(1-delta0);
}

} // end loop over cases


}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

vector[3] thetaalt;	
vector[3] thetaalts;
real pialt;

real inconsistj[M];	
real inconsistjj[M];	
real inconsist;	
real pt1;
real pt2;

theta = inv_logit(logittheta);

for (j in 1:M){

llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);

// simulate probability of different outcome under alternative panels
        for (jj in 1:M){ 
            thetaalt[1] = inv_logit(logittheta[Ij[jj,1]] + phi[Lj[j]]);
            thetaalt[2] = inv_logit(logittheta[Ij[jj,2]] + phi[Lj[j]]);
            thetaalt[3] = inv_logit(logittheta[Ij[jj,3]] + phi[Lj[j]]);
            thetaalts = sort_asc(thetaalt);
            // Define mixture model from submodels, split by covariate
            if (SC[j]==1) { // if cov = 1
            pialt = thetaalt[1]*delta1 + thetaalts[2]*(1-delta1);
            } else {
            pialt = thetaalt[1]*delta0 + thetaalts[2]*(1-delta0);
            }
            
            pt1 = (Y[j]*(pi[j] > pialt)*((pi[j] - pialt)/(pi[j])));
            pt2 = ((1-Y[j])*(pi[j] < pialt)*((pialt - pi[j])/(1 - pi[j])));
            inconsistjj[jj] = (pt1 + pt2);
            
        }
    inconsistj[j] = mean(inconsistjj);
    }
ll = sum(llj);
inconsist = mean(inconsistj);


}

'



stancode_covariate_mixture_independent = '

data {
int<lower=1> M;	
int<lower=1> N;
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
real<lower=0,upper=1> delta;
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector<lower=0,upper=1>[2] pimat[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M){	// loop over cases

// Define chair submodel
pimat[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);

// Define median submodel
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);
thetamats3[j] = rep_vector(0.0,3);        
pimat[j,2] = 0.0;

thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);

// probability of at least two yes votes, assuming judges make conditionally independent evaluations
pimat[j,2] = thetamat3[j,1]*thetamat3[j,2]*(1-thetamat3[j,3]) + thetamat3[j,1]*(1-thetamat3[j,2])*thetamat3[j,3] + (1-thetamat3[j,1])*thetamat3[j,2]*thetamat3[j,3] + thetamat3[j,1]*thetamat3[j,2]*thetamat3[j,3]; 	

// Define mixture model from submodels

pi[j] = pimat[j,1]*delta + pimat[j,2]*(1-delta);


} // end loop over cases


}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;


vector[3] thetaalt;	
real pialt;
real inconsistj[M];	
real inconsistjj[M];	
real inconsist;	

theta = inv_logit(logittheta);

for (j in 1:M){

llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);


// simulate probability of different outcome under alternative panels
        for (jj in 1:M){ 
            thetaalt[1] = inv_logit(logittheta[Ij[jj,1]] + phi[Lj[j]]);
            thetaalt[2] = inv_logit(logittheta[Ij[jj,2]] + phi[Lj[j]]);
            thetaalt[3] = inv_logit(logittheta[Ij[jj,3]] + phi[Lj[j]]);
            
            // probability of at least two yes votes, assuming judges make conditionally independent evaluations
            pialt = thetaalt[1]*thetaalt[2]*(1-thetaalt[3]) + thetaalt[1]*(1-thetaalt[2])*thetaalt[3] + (1-thetaalt[1])*thetaalt[2]*thetaalt[3] + thetaalt[1]*thetaalt[2]*thetaalt[3]; 	
            pialt = thetaalt[1]*delta + pialt*(1-delta); // define predicted probability from sub-models
            inconsistjj[jj] = (Y[j]*(1-pialt) + (1-Y[j])*pialt);
            
        }
    inconsistj[j] = mean(inconsistjj);
}
ll = sum(llj);	
inconsist = mean(inconsistj);


}

'


###################
## median models ##
###################


stancode_covariate_median = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	


for (j in 1:M){	// loop over cases

// Define median submodel
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);
thetamats3[j] = rep_vector(0.0,3);        
pi[j] = 0.0;


// if all judges are known
thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);
thetamats3[j] = sort_asc(thetamat3[j]);
pi[j] = thetamats3[j,2]; 		

} // end loop over cases

}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

vector[3] thetaalt;	
vector[3] thetaalts;
real pialt;
real inconsistj[M];	
real inconsistjj[M];	
real inconsist;	
real pt1;
real pt2;

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);


// simulate probability of different outcome under alternative panels
        for (jj in 1:M){ 
            thetaalt[1] = inv_logit(logittheta[Ij[jj,1]] + phi[Lj[j]]);
            thetaalt[2] = inv_logit(logittheta[Ij[jj,2]] + phi[Lj[j]]);
            thetaalt[3] = inv_logit(logittheta[Ij[jj,3]] + phi[Lj[j]]);
            thetaalts = sort_asc(thetaalt);
            pialt = thetaalts[2];
            pt1 = (Y[j]*(pi[j] > pialt)*((pi[j] - pialt)/(pi[j])));
            pt2 = ((1-Y[j])*(pi[j] < pialt)*((pialt - pi[j])/(1 - pi[j])));
            inconsistjj[jj] = (pt1 + pt2);
            
        }
    inconsistj[j] = mean(inconsistjj);
}
ll = sum(llj);
inconsist = mean(inconsistj);

}
'





stancode_covariate_median_independent = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];	


for (j in 1:M){	// loop over cases

// Define median submodel
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);     
pi[j] = 0.0;

thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);

// probability of at least two yes votes, assuming judges make conditionally independent evaluations
pi[j] = thetamat3[j,1]*thetamat3[j,2]*(1-thetamat3[j,3]) + thetamat3[j,1]*(1-thetamat3[j,2])*thetamat3[j,3] + (1-thetamat3[j,1])*thetamat3[j,2]*thetamat3[j,3] + thetamat3[j,1]*thetamat3[j,2]*thetamat3[j,3]; 		


} // end loop over cases

}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

vector[3] thetaalt;	
vector[3] thetaalts;
real pialt;
real inconsistj[M];	
real inconsistjj[M];	
real inconsist;	


theta = inv_logit(logittheta);

for (j in 1:M){

llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);


// simulate probability of different outcome under alternative panels
  for (jj in 1:M){ 
    thetaalt[1] = inv_logit(logittheta[Ij[jj,1]] + phi[Lj[j]]);
    thetaalt[2] = inv_logit(logittheta[Ij[jj,2]] + phi[Lj[j]]);
    thetaalt[3] = inv_logit(logittheta[Ij[jj,3]] + phi[Lj[j]]);

    // probability of at least two yes votes, assuming judges make conditionally independent evaluations
    pialt = thetaalt[1]*thetaalt[2]*(1-thetaalt[3]) + thetaalt[1]*(1-thetaalt[2])*thetaalt[3] + (1-thetaalt[1])*thetaalt[2]*thetaalt[3] + thetaalt[1]*thetaalt[2]*thetaalt[3]; 		

    inconsistjj[jj] = (Y[j]*(1-pialt) + (1-Y[j])*pialt);
    
  }
  inconsistj[j] = mean(inconsistjj);
}
ll = sum(llj);	
inconsist = mean(inconsistj);

}
'



###############
## max model ##
###############


stancode_covariate_max = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M){	// loop over cases

// Define max submodel		
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);
thetamats3[j] = rep_vector(0.0,3);        
pi[j] = 0.0;

// if all judges are known
thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);
thetamats3[j] = sort_asc(thetamat3[j]);
pi[j] = thetamats3[j,3];

} // end loop over cases

}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

theta = inv_logit(logittheta);

for (j in 1:M){
    llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);
}
ll = sum(llj);
}
'

###############
## min model ##
###############



stancode_covariate_min = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M){	// loop over cases

// Define min submodel		
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);
thetamats3[j] = rep_vector(0.0,3);        
pi[j] = 0.0;

// if all judges are known
thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);
thetamats3[j] = sort_asc(thetamat3[j]);
pi[j] = thetamats3[j,1];


} // end loop over cases


}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

theta = inv_logit(logittheta);

for (j in 1:M){
    llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);
}
ll = sum(llj);
}
'



##################
## chair models ##
##################


stancode_covariate_chair = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M){	// loop over cases

// Define chair submodel
pi[j] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);

} // end loop over cases



}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

real pialt;
real inconsistj[M];	
real inconsistjj[M];	
real inconsist;	
real pt1;
real pt2;

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);

// simulate probability of different outcome under alternative panels
        for (jj in 1:M){ 
            pialt = inv_logit(logittheta[Ij[jj,1]] + phi[Lj[j]]);
            pt1 = (Y[j]*(pi[j] > pialt)*((pi[j] - pialt)/(pi[j])));
            pt2 = ((1-Y[j])*(pi[j] < pialt)*((pialt - pi[j])/(1 - pi[j])));
            inconsistjj[jj] = (pt1 + pt2);
            
        }
    inconsistj[j] = mean(inconsistjj);
}
ll = sum(llj);
inconsist = mean(inconsistj);

}
'


stancode_covariate_chair_independent = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M){	// loop over cases

// Define chair submodel
pi[j] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);

} // end loop over cases



}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

real pialt;
real inconsistj[M];	
real inconsistjj[M];	
real inconsist;	

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);


// simulate probability of different outcome under alternative panels
        for (jj in 1:M){ 
            pialt = inv_logit(logittheta[Ij[jj,1]] + phi[Lj[j]]);
           inconsistjj[jj] = (Y[j]*(1-pialt) + (1-Y[j])*pialt);
            
        }
    inconsistj[j] = mean(inconsistjj);
}
ll = sum(llj);	
inconsist = mean(inconsistj);

}
'




##################
## second model ##
##################


stancode_covariate_second = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M) pi[j] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);

}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]); 
}
ll = sum(llj);

}
'

#################
## third model ##
#################



stancode_covariate_third = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	

for (j in 1:M) pi[j] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);

}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);
}
ll = sum(llj);

}
'

################
## null model ##
################


stancode_covariate_null = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real<lower=0> sigmaphi;
real logittheta;
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];

for (j in 1:M) pi[j] = inv_logit(logittheta + phi[Lj[j]]);

}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta;
real llj[M];	
real ll;

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]); 
}
ll = sum(llj);

}
'


################
## mean model ##
################

stancode_covariate_mean = '

data {
int<lower=1> M;	
int<lower=1> N;	
int<lower=1> K;	
int<lower=0,upper=1> Y[M];		
int<lower=0,upper=N> Ij[M,3]; 
int<lower=1,upper=K> Lj[M];	
}

parameters { 
real mutheta;
real<lower=0> sigmatheta;
real<lower=0> sigmaphi;
real logittheta[N];
vector[K] phi;
}

transformed parameters { 

real<lower=0,upper=1> pi[M];
vector[3] thetamat3[M];
vector[3] thetamats3[M];	


for (j in 1:M){	// loop over cases

// Define median submodel
// initialize objects to zero
thetamat3[j] = rep_vector(0.0,3);
thetamats3[j] = rep_vector(0.0,3);        
pi[j] = 0.0;

// if all judges are known
thetamat3[j,1] = inv_logit(logittheta[Ij[j,1]] + phi[Lj[j]]);
thetamat3[j,2] = inv_logit(logittheta[Ij[j,2]] + phi[Lj[j]]);
thetamat3[j,3] = inv_logit(logittheta[Ij[j,3]] + phi[Lj[j]]);
pi[j] = mean(thetamat3[j]); 		// mean model


} // end loop over cases

}	

model {		  	

for (j in 1:M) Y[j] ~ bernoulli(pi[j]);

for (i in 1:N) logittheta[i] ~ normal(mutheta,sigmatheta);

for (k in 1:K) phi[k] ~ normal(0,sigmaphi);

}

generated quantities {

real<lower=0,upper=1> theta[N];
real llj[M];	
real ll;

theta = inv_logit(logittheta);

for (j in 1:M){
llj[j] = Y[j]*log(pi[j]) + (1-Y[j])*log(1-pi[j]);
}
ll = sum(llj);

}
'








